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Abstract 

In this letter we address issues involved in quarkonium production near the boundaries of 
phase space. It is shown that higher-order non-perturbative contributions are enhanced 
in this kinematic region and lead to a breakdown of the non-relativistic (NRQCD) ex- 
pansion. This breakdown is a consequence of sensitivity to the kinematics of soft gluon 
radiation and to the difference between partonic and hadronic phase space. We show 
how these large corrections can be resummed giving the dominant contribution to the 
cross section. The resummation leads to the introduction of non-perturbative, universal 
distribution functions. We discuss the importance of these shape functions for several 
observables, in particular the energy distribution of photo-produced J/ip close to the 
endpoint. 



Non-relativistic effective field theory (NRQCD) M for quarkonium decays and pro- 
duction has shed new light on many aspects of the experimental data. In particular, 
the inclusion of color-octet decay and production channels rectified some of the main 
theoretical and phenomenological deficiencies of the color-singlet model and could ex- 
plain the so-called '^'-anomaly' in pp collisions at the Tevatron || . Despite this success, 
the verification of the universality of long-distance parameters in NRQCD through the 
comparison of different production processes turns out to be difficult due to significant 
theoretical uncertainties of various kinds. A large uncertainty stems from the treatment 
of the kinematics of the hadronization process, which is especially important near the 
boundary of phase space where keeping the leading order terms in the NRQCD expansion 
is not valid. Here we concern ourselves exactly with this problem. A partial resummation 
of the NRQCD expansion leads us to introduce universal shape functions that parame- 
terize quarkonium formation at its kinematic limits. We discuss the effect of these shape 
functions on quarkonium production rates in a variety of production processes. 

In the NRQCD approach the inclusive quarkonium production process factorizes into 
the production of a (heavy) quark-anti-quark pair at small relative momentum 

initial state — > QQ[n] + X, (1) 
followed by the 'hadronization' of the QQ pair into a quarkonium H and light hadrons 

QQ[n]^H + X. (2) 
The cross section may be written as 

dcr H = H da QQ[n]+x(°n ), (3) 

n 

where g^qq^j+x is the short distance cross section to produce a QQ state labeled by n 
plus X, and the general form of the NRQCD matrix element (O^ ) is given by 

{OS) = E (Q\^nX\H + X)(H + X| X t r >|0>. (4) 
x 

T n (T' n ) is a matrix in color and spin and may contain derivatives and, in principle, gluon 
and light quark fields. 

The short-distance cross section is insensitive to quarkonium binding by construc- 
tion and therefore can depend only on the heavy quark mass and kinematic invariants 
constructed from parton momenta. The quarkonium mass never enters explicitly, since 
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any non-perturbative effect is parametrized by one of the matrix elements (O^ )■ As a 
further consequence of factorization, the momentum taken by the light hadronic final 
state X in (0) is neglected at lowest order in the velocity expansion [Ij. This is consistent 
since the energies involved in (|2]) are small, of order rriQV 2 , in the QQ rest frame, and 
because the process is inclusive over the light hadronic final state. Thus, given that the 
difference between partonic and hadronic kinematics is higher order in v 2 , shifting say 
a partonic to a hadronic endpoint, does not in general improve the accuracy of a calcu- 
lation if other effects of the same order in the velocity expansion are left out. However, 
there are circumstances, namely near the boundaries of phase space, where effects higher 
order in v 2 related to kinematics are enhanced and lead to a breakdown of the velocity 
expansion. In these cases it becomes necessary to include these effects in a systematic 
fashion. The breakdown of the velocity expansion is closely related to the fact that near 
the boundaries of phase space the details of the hadronization process (0) are probed 
at a scale smaller than rrigv 2 . This happens, for instance, for quarkonium production 
close to threshold, s —>■ M^, or in J/ip photo-production close to the point where all the 
photon's energy is transferred to the J/ip (z = Ej/^/E^ — > 1 in the proton rest frame). 
In these cases, to make the NRQCD expansion convergent, one must smear sufficiently 
(in s and z, respectively) so as not to probe the details of hadronization. However, it is 
possible to reduce the necessary smearing region by resumming a class of leading twist 
operators into universal shape functions, thus extending the range of applicability of 
NRQCD. Resummations of related nature have been introduced to treat energy spectra 
in semileptonic or radiative B meson || and quarkonium Q decays. 

Each kinematic situation must be dealt with slightly differently, although the general 
idea is the same. Consider some general production process A + B — > H + X. To 
determine the matching coefficients ^o"QQ[ n ] + x at lowest order we calculate the inclusive 
rate in full QCD to produce a free quark anti-quark pair in a state n. The momenta of 
the quark and anti-quark are defined as 

p„ = *f + (Ah), p, = y + W/" ( 5 ) 

respectively. A is the Lorentz transformation which boosts from the frame where 
is given by (2m<g, 0) to whatever frame one chooses for the calculation. In the P-rest 
frame k = ( \Jff + rriQ — mQ,li) so that p 2 = p 2 = nig. The production cross section 
is then expanded in the small momenta U ~ 0(rriQv), and powers of momentum are 
identified with derivatives acting on heavy quark fields in NRQCD. Thus, a factor of 
relative momentum In — In corresponds to (ipHDix) and a factor of center-of-mass (cms) 
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momentum (of the QQ pair in the P-rest frame) lu+hi is identified with a total derivative 
iDify'x)- Four-fermion operators with total derivatives on fermion bilinears are usually 
ignored since they are of higher order in the non-relativistic expansion than the relative 
momentum operators. They were first written down in || and were later shown to be 
of dynamical importance in radiative quarkonium decay near the endpoint |4], |J. (Note 
that using the equations of motion Z10 + /20 can be identified with the total time derivative 

It is precisely these cms-derivative operators that cause the leading singular contri- 
butions in the NRQCD expansion in a general situation of constrained phase space. The 
phase space measure for the process, initial state — > Q{p) + Q{p) + X, is given by 

dps = , d "S TT , , 5*(yt l -Yk i -(p + A(h + 1 2 ))) , (6) 

(2tt) 3 2P i (27i) 3 2k V « / 

where the t\ are the incoming momenta and ki are the momenta of final state particles 
other than the QQ pair. The matching strategy then tells us to expand this delta- 
function, as well as the amplitude squared, in powers of the momenta li and l^. The delta- 
function depends explicitly on the cms momentum l\ + l 2 and implicitly on the relative 
momentum through the zero-components of ij. The expansion of the delta-function 
results in a series in v 2 , but with increasingly singular distributions at the boundary 
of partonic phase space defined by the delta-function with li set to zero. Normally, 
the phase space is integrated with a smooth function, and the terms coming from the 
expansion of the delta function converge rapidly, provided v 2 is small. However, if the 
weighting function has support mainly in a region of order erriQ near the boundary of 
phase space, the expansion parameter is v 2 /e. Thus e ^> v 2 is required for convergence. 
On the other hand, resumming the leading terms (v 2 /e) n to all orders allows us to take 
e ~ v 2 . Provided that the squared amplitude for the production process is non-singular 
at the boundary of phase space, the resummation of the leading singular contributions 
is easily accomplished. Note that even after resummation we can not consider e<Cti 2 , 
because all subleading terms v 2m {v 2 / e) n become important as well in this region. Since 
the details depend on the particular kinematic situation, we shall show how this works 
in several processes of topical interest. 

Let us begin by considering the case of gluon fragmentation into ip {J/4> or 4>') 
through a QQ pair in a 3 S\ color-octet state, which is now regarded as the dominant 
source of direct ip production at large transverse momentum in hadron collisions . The 
octet QQ pair evolves non-perturbatively into the ip, a process necessarily accompanied 
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by soft gluon radiation. At lowest order in v 2 the momentum of the soft gluon radiation 
can be neglected and and the fragmentation function is given by 

The hat reminds us that the variable z = P + /k + is defined in terms of the light-cone 
+- components of parton momenta, that is, P + refers to P in ([5]) rather than the ip- 
momentum. k is the momentum of the fragmenting gluon. As noted in Jj]] this fragmen- 
tation function is folded with a gluon production cross section, which, for a given value 
of the ip transverse momentum, is a rapidly varying function of z. To be precise 

da r dz 



d j ~ K{z,p t )D^ /g {z, Pt ), (8) 

where K(z,p t ) ~ z 5 in the region p t ~ 15 GeV. Consequently, the integral (H) depends 
sensitively on the shape of the fragmentation function near z = 1 and neglecting the kine- 
matics of soft gluon emission is not a good approximation. Intuitively, one expects soft 
gluon radiation to soften the fragmentation function^] (|7|) and therefore to reduce da/dp t 
M. However, this intuition relies on interpreting z as the quarkonium +- momentum 
fraction, different from the definition of z given above. 

Now the existence of cms momentum I = lx + l 2 in the matching relations can be 
interpreted in the quarkonium rest frame as recoil of the QQ pair against the emitted 
soft gluons. Indeed, since the matrix element for g* — > QQ is evidently non-singular at 
z = 1, the leading singular contributions to the coefficient functions of higher- dimension 
operators in the NRQCD expansion simply follow from (||). After integrating over phase 
space except for the +- direction, this leads to the shift 

6{l-z)^6(l-z-A + ^/P + ). (9) 

In terms of operators in the NRQCD expansion, expanding the modified delta-function 
in A +fJi l^/P + = /+/ (2mg), we obtain the series 

5>|^T^ + + X\(in ■ D) m (xV.T^) |0), (10) 
x 



1 This softening occurs before the onset of Altarelli-Parisi evolution that further softens the fragmen- 
tation function at the scale pt- 
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where n is the light like vector (1,0,0,-1), D = D/(2m c ) and 5^ m '(l — z) denotes 
the m'th derivative of the delta-function. The first term in the sum reproduces ([?]). 
Introducing the shape function 

F4 3 S{ 8 % + ) = Y.^°iT A x\4, + X)^ + X\5(y + -in-D) ( X V,:rty) |0), (11) 
x 

we rewrite (|I0|) as 

D^ g (z,2m c ) = J dy + 5(l-z-y + )F^s[ 8) ](y + ). (12) 

Note that F i ,[ 3 S{ 8) ](y + ) is a non-perturbative distribution function that can not be calcu- 
lated except within models. But the factorization as expressed by (|ll|) and ([12]) implies 
that such shape functions are process- independent. Moreover, F^S^^y-^) has support 
mainly for \y + \ < v 2 . We assume ttiqV 3> vtiqv 2 ~ Aqcd for NRQCD power counting, in 
which case a cms-derivative acting on a quark bilinear scales as m,QV 2 , the momentum of 
dynamical gluons in NRQCD.0 It is suggestive to interpret F^p[ 3 S^}(y + ) as the probabil- 
ity for soft gluons to carry off a light-cone momentum fraction y + in the hadronization of 
the color-octet QQ pair. Although useful, this interpretation is not quite exact, because, 
as will be seen in more detail below, such distribution functions are non-trivial even 
for color-singlet QQ states, when no soft gluon emission is required. In such cases, the 
shape functions account for the difference between quark and quarkonium masses and 
the corresponding difference in the definition of kinematic variables. 

The ip production rate at fixed p t follows from integrating the fragmentation function 
with a function that roughly varies as i 4 , see (|8]). Eq. ([T0| ) implies that the NRQCD 
expansion parameter is then 4t> 2 ~ 1 rather than v 2 . Inclusion of the effects due to the 
shape function F 1 p[ 3 s[ 8 ' > }(y + ) can therefore alter the ip production cross section of order 
unity, an effect that should be kept in mind when (Og ( 3 5*i)) extracted from Tevatron 
data is compared with the extraction from another production process. Unfortunately, 
without a non-perturbative calculation of the shape function, no detailed prediction is 
possible. We note that the prediction that the ip is transversely polarized at large pt [JTO] 



in hadron collisions remains unaffected by the resummation procedure described here. 

(For perturbative corrections to this prediction see [p|] .) 

2 In the construction of ]8| this can be seen from the fact that the matrix element with n cms 
derivatives vanishes unless one inserts an n'th multipolc operator resulting in a an additional 1/c™ 
factor above those coming from naive dimensional analysis. 
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Our second example concerns the total hadro-production cross section of a quarko- 
nium H at fixed target energies. In this case we write the cross section as 

1 

<Jh = J2 J dx idx 2 fi/A(x 1 )f j/B (x 2 ) o(ij ->• H + X) , (13) 

i,j o 

The sum extends over all partons in the colliding hadrons and fa/ a, etc. denote the 
corresponding parton distribution functions. The parton cross section a is a distribution 
that is weighted with the parton densities. Let us consider the leading parton process 
i + j — > Q(p) + Q(p)- Since its matrix element is again non-singular at the boundaries 
of phase space, it is sufficient to retain in (^j) and to set k = in the matrix element, 
in order to obtain the most singular coefficient functions to all orders in v 2 . The parton 
cross section may then be written as 



7r 3 a 2 



(2m Q ) 



£(0|^r n x|# + X)(H + X\5( Xl x 2 s - 4m 2 Q ~ 2 ( k i + k i) " A/ ) xKm, (14) 



x 



where kij are the four-momenta of the incoming partons, s the hadron center-of-mass 
energy and I = iD. The sum over n extends over all possible operators and is truncated 
according to the non-relativistic power counting rules. For the two most important 
production channels in S-wave quarkonium production C^S®} = 5/12, C«M*>] = 
35/(12m£) 0, 0. Using 2{k, + k,) ■ Al » 2P ■ Al = 4m Q l° and expanding the delta- 



function gives the series of higher-dimension operators 



[2m Q fs V m ml s 



in 



U^r nX \H + x)(H + x\[ 8 -^^\ (xKri\o), (15) 



X 



which accounts for all most singular contributions at x±x 2 = 4m 2 -)/ s. As previously, 
D = D/(2mo). Note that as opposed to flTU|), the higher- dimension operators involve 
time rather than +- derivatives. Accordingly we define a shape function in energy fraction 
Ve 

F H [n}(y E ) = J2(^ T nX\H + X)(H + X\6(y E - W ) (xK^) |0>- (16) 



x 



We emphasize that Fh\p\{ve) and Fjj[n](y+) in (|TT1) are different non-perturbative dis- 



tributions. Up to the caveat mentioned above, F H [n](y E ) can be interpreted, in the 
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quarkonium rest frame, as the distribution of energy fraction taken by soft gluons during 
the transition from the short- distance heavy quark pair in a state n into the quarkonium 
H. 

Consider now the case where the QQ pair is in the same color-singlet state as the 
quarkonium H. Up to corrections suppressed in v 2 , we can use the vacuum saturation 
approximation and drop the sum over X in (fL6|). The derivative iD Q produces a factor 
of the binding energy and we obtain 

F H [n](y E ) = 5{y E - (Q\^T nX \H) (H\^T> n m , (17) 

where M# is the quarkonium mass. Consequently, the ^-integral that enters the pro- 
duction cross section (see ( p0[ ) below) can be done: 

Jdy E F H [n](y E ) 5(xix 2 s - 4m 2 Q - 8m 2 Q y E ) = 5{x l x 2 s - Am 2 Q - Am Q (M H - 2m Q )) 

■(o\^r nX \H)(H\ x Km- (18) 

Since NL\ = Am 2 Q+AmQ(M H —2mQ)+0{y A ), we see that the sole effect of the distribution 
function in the vacuum saturation approximation is to shift the unphysical partonic 
boundary of phase space to the hadronic one. The inclusion of sub-leading terms would 
complete the transmutation 

5{xi% 2 s -4mg) -y 5(x!X 2 s - M 2 H ). (19) 

As we previously mentioned, performing this resummation does not in general im- 
prove the accuracy of the calculation, given that other contributions of the same or- 
der in v 2 have been left out. However, if the corrections from the expansion of the 
delta-function are enhanced because the dominant contribution to the cross section 
comes from a region close to the boundary of phase space, then we have indeed im- 
proved the situation. Introducing the variable r = x\x 2 and the parton luminosities 
Lij(r) = f^: dx/x fi/A(x)fj/B(i~/x), the hadro-production cross section can be expressed 
as 

° H = ?3 q) £ /W)E^H Jdy E F H [n](y E )5(TS-4m 2 Q -8m 2 Q y E ). 



h3 



(20) 

At high energies s ^> 4mg, the gluon-gluon fusion channel dominates and we see that 
the partonic cross section is weighted by a gluon luminosity, that, because of the small-x 
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behaviour of the gluon distribution, increases rapidly as r decreases. Thus, a systematic 
shift in t due to an asymmetric shape function Fh\p\{ve) (such as the example (|17|) 
above) can affect the magnitude of the total cross section in an important way. Indeed, 
one finds in the leading order in v 2 predictions that the difference between quark and 
quarkonium masses in the phase space delta-function introduces a normalization uncer- 



tainty of up to a factor two for ip' and x production. As noted in [|12| , the effect should 
be even more pronounced in color-octet mechanisms as, due to soft gluon radiation, the 
final state invariant mass is always larger than M\. 

As our final, and perhaps most interesting, example we examine photo-production of 
J /ip, in particular the ^-distribution, where z = p-Pj/^/p-k^ and p is the proton momen- 
tum. 3 (z = Ej/^/Ery in the proton rest frame.) To exclude higher-twist contributions 
from diffractive J ftp production, which can not be accommodated in the leading-twist 
NRQCD factorization formalism, it is preferable to consider inelastic J/ip production 
and to impose z < 0.9 and P t > 1 GeV. Inelastic J /ip photo-production is often per- 
ceived as a problem for the NRQCD theory of onium production, because the color-octet 
contributions to this process rise rapidly with z ||14|| , which is not supported by the data. 
On the other hand octet production is accompanied by gluon radiation, which carries 
away an energy fraction of order v 2 and makes it highly unlikely to reach a value of z 
close to one [T^]. Technically, the necessary smearing of the energy distribution is again 
achieved through a shape function, as we shall see, although the situation is slightly 
more complicated for the 2 — > 3 parton processes than for the 2^2 parton processes 
considered before. 

Before addressing this issue, let us briefly comment on the 2 — > 2 parton processes in 
photo-production, which are kinematically restricted to z — p ■ P/p ■ k 7 = 1 and P t = 0. 
We emphasize again that the calculation in NRQCD refers to partonic variables like z 
with P defined in (Q) rather than z. If we are interested in the total cross section then 
the analysis goes through almost exactly as in the case of hadro-production. In fact the 
same shape function (|16D enters. However, the total cross section is of marginal interest 
experimentally. More interesting are the differential distributions. For the gluon fusion 
process 



l 

/ \ f , , ,9 \ ^ ^ynr n ^ r n / 

Pt) 



3 Because we will be concerned with the large-z region, we do not consider resolved photon interactions 
in the following. 
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6(sx - 4mJ ~ 8m QVE) 6(1- z- y+) 6®{P t - p t ), (21) 

where 

F H [n}(y E ,y + ,p t ) = Y,(0\^nX\H + X) 



x 



(H + X\6(y + - in ■ D)5(y E - iD»)5{p T - iD T ) {^Vj)) |0>. (22) 

Such a mult i- dimensional distribution is much too complicated to be useful in practice, 
and we do not pursue it further here. Qualitatively, this distribution causes smearing in 
z and transverse momentum over a region (for Q = c) 

5z^v 2 « 0.25 - 0.3, 5P t w m c v 2 w 0.5 GeV. (23) 

Unless the endpoint region is averaged over an interval of this size, the NRQCD expansion 
fails to be convergent, even if the above distribution function is taken into account. 
Smearing over a range much larger than (|23D , makes the corrections considered here 
irrelevant and one can return to the ordinary velocity expansion. Note that because of 
the sizeable smearing in z, part of the 2 — > 2 parton contributions that are normally 
considered to be localized at z = 1 actually leaks into the inelastic region z < 0.9. This 
leakage can be eliminated by a transverse momentum cut P t ^> m c v 2 . Note also that the 
smearing in the transverse direction is of the same order of magnitude as the smearing 
caused by intrinsic transverse momentum of the gluon in the proton. 

Let us return now to inelastic photo-production of quarkonium and consider 7 + g —>■ 
Qip) + Qip) + g(k). It is straightforward to obtain 

da 



■£= J dP 2 JdxS(x,z,p 2 )Y,(o\^r nX \H + x) 

t->2 n X 



(/If 2(1 _ £\ _|_ p2 \ 
s(l - z)x 1 1 - 2k ■ AtDj (xKiP) |0>, (24) 

S(x } z,P 2 ) = —^—f g/p ( x )\M n \ 2 , 
lbirzxs 

where s is the cms energy, M = 2rriQ, and M. n is the matrix element to produce the 
QQ pair in a state n. If we neglect the covariant derivative D in the delta-function, the 
standard kinematic relations for the inelastic processes considered in Jbl are recovered. 



To proceed, we have to extract the dependence of 2k ■ AD on z and Pt- To this end we 
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note that 



J2(0\^ nX \H + X)(H + XKiD^) . . . (zD,J (xK^) |0> 



x 



An fa] w^... w Mm + g mix -terms (25) 



in a frame where the QQ-pa.ii moves with four-velocity w. Since k is a light-like four- 
vector, terms involving the metric tensor <? MiA i do not contribute and we obtain, in the 
rest frame of P (see (|5|)), 

J2(0\^T nX \H + X}(H + X\(2k-MDr(xK^) |0> = A»[n] (2fc^) m . (26) 

Using 

we can write Q23D as 



2P-k _ M 2 (l-z) 2 + P 2 
Kest ~ M ~ Mz(l-z) ' K ' 



-£= J dP t 2 JdxJ dy + S{x, z, P, 2 ) F H [n]{y, 



^(l-^- ^ 1 -/)^ - ^-!^^ ), (28) 

where the shape function i 7 #[n](2/+) is defined as in (p] ), generalized to an arbitrary 
QQ state n. Now let us examine the delta-function. For any non-zero Pt we see that 
the expansion parameter close to z — 1 is y+/(l — z) ~ t> 2 /(l — 5). Thus, the NRQCD 
expansion breaks down for 1 — z ~ 0(v 2 ), because higher-order terms in v 2 grow more 
and more rapidly as z — > 1.Q Consequently, the NRQCD factorization approach makes no 
prediction in the endpoint region and the discrepancy between leading order predictions 
|14| and data in this region does not allow us to draw any conclusion on the relevance 



of color-octet contributions to photo-production. If one averages the z-distribution over 
a sufficiently large region containing the endpoint, the octet mechanisms contribute 
significantly to this average. However, the characteristic shape information is then lost 
and one has to deal with the more difficult and uncertain question of whether the absolute 



4 The integral over z still exists, because for any fixed Pt, z can never reach one. However, 1 — z max = 
0(M 2 /s, P t 2 / s) is very small in a high energy collision. We also mention that in the case considered 
here, the expansion of the partonic amplitude in l\ and li can also produce terms more singular in 
1/(1 — z) in higher orders in li. 
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magnitude of the cross section requires the presence of octet contributions and whether 
their magnitude is consistent with other production processes. Almost identical remarks 
apply to the endpoint region in the J /ip energy distribution in the process B — > J/ip + X 
or e + e~ — > J/if) + X. 

In conclusion, we have shown how the NRQCD expansion of quarkonium production 
processes breaks down in regions of phase space, where the prediction is sensitive either 
to the difference between parton and hadron kinematics or to momentum carried away 
by light hadrons in the hadronization process (0). The amount of smearing which is ne- 
cessitated by this breakdown can be reduced by introducing universal shape functions in 
NRQCD as illustrated in this letter. Being universal, these functions could in principle 
be extracted from one production process and used for another. However, the existence 
of such functions for every QQ state n at short distances, together with functions in 
different variables (such as y + and de) makes this difficult in practice. At this stage 
it seems worthwhile to use the operator expressions for the shape functions to invent 
phenomenological models for these functions, which are consistent with the properties 
of NRQCD. 
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